clear
rng default

load Cw_R1
C_21_w_R1=C_21_w;
C_32_w_R1=C_32_w;
C_43_w_R1=C_43_w;

load Cw_R2
C_21_w_R2=C_21_w;
C_32_w_R2=C_32_w;
C_43_w_R2=C_43_w;


all_years=1940:1:1967;
n_years=size(all_years,2);
ntypes=4;

C_21_w_Gumbel=zeros(n_years,1);
C_32_w_Gumbel=zeros(n_years,1);
C_43_w_Gumbel=zeros(n_years,1);
C_21_w_Gumbel_original=zeros(n_years,1);
C_32_w_Gumbel_original=zeros(n_years,1);
C_43_w_Gumbel_original=zeros(n_years,1);



for j=1:n_years
    
year=all_years(j);

% Probability distribution of Ychosen|X
condprobsmen_temp=load('probs_allyears_data.mat',['condprobsmen' num2str(year) '_140']);
condprobsmen=condprobsmen_temp.(['condprobsmen' num2str(year) '_140']);
PYX_cond=[condprobsmen(:,2:end) condprobsmen(:,1)]; %4x5
%Note: PYX_cond(x,y)=Probability to marry woman y born in any year conditional on being man x born in 1940

% Probability distribution of Xchosen|Y
condprobswomen_temp=load('probs_allyears_data.mat',['condprobswomen' num2str(year+1) '_140']);
condprobswomen=condprobswomen_temp.(['condprobswomen' num2str(year+1) '_140']);
PXY_cond=[condprobswomen(:,2:end) condprobswomen(:,1)];  %4x5
%Note: PXY_cond(y,x)=Probability to marry man x born in any year conditional on being woman y born in 1941

% Probability distribution of X
probsmen_temp=load('probs_allyears_data.mat',['probsmen' num2str(year) '_140']);
PX=probsmen_temp.(['probsmen' num2str(year) '_140']);
%PX(x)=Probability of being man x conditional on being born in 1940

% Probability distribution of Y
probswomen_temp=load('probs_allyears_data.mat',['probswomen' num2str(year+1) '_140']);
PY=probswomen_temp.(['probswomen' num2str(year+1) '_140']);
%PY(y)=Probability of being woman y conditional on being born in 1941

[~, V_Gumbel, ~, ~,~,~]=param_Gumbel_closedform(PYX_cond, PXY_cond, ntypes, ntypes);


Wmarital_systematic=zeros(1,ntypes-1);
for h=1:ntypes-1
    Wmarital_systematic(j,h)=sum(PXY_cond(h+1,1:end-1).*(V_Gumbel(:,h+1).'),2)-sum(PXY_cond(h,1:end-1).*(V_Gumbel(:,h).'),2);
end

Wmarital=zeros(1,ntypes-1);
for h=1:ntypes-1
    Wmarital(j,h)=(log(exp(V_Gumbel(1,h+1))+exp(V_Gumbel(2,h+1))+exp(V_Gumbel(3,h+1))+exp(V_Gumbel(4,h+1))+1)+ 0.5772)- ...
                       (log(exp(V_Gumbel(1,h))+exp(V_Gumbel(2,h))+exp(V_Gumbel(3,h))+exp(V_Gumbel(4,h))+1)+ 0.5772);
end


C_21_w_Gumbel(j)=Wmarital_systematic(j,1);
C_32_w_Gumbel(j)=Wmarital_systematic(j,2);
C_43_w_Gumbel(j)=Wmarital_systematic(j,3);

C_21_w_Gumbel_original(j)=Wmarital(j,1);
C_32_w_Gumbel_original(j)=Wmarital(j,2);
C_43_w_Gumbel_original(j)=Wmarital(j,3);

end




%% Figure 5 in the paper

figure
plot(all_years, C_21_w_Gumbel, 'b', 'LineWidth',3)
hold on
plot(all_years, C_21_w_Gumbel_original, 'k', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
box on
xlim([1940 1967])
ylim([ -1 1 ])
set(gca,'ytick', -1:0.5:1 , 'FontSize',15)
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{21}(V) and \\Delta_{21}(V)'),'FontSize',15)
saveas(gcf, 'C_21_w_original.jpg')


figure
plot(all_years, C_32_w_Gumbel, 'b', 'LineWidth',3)
hold on
plot(all_years, C_32_w_Gumbel_original, 'k', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
box on
xlim([1940 1967])
ylim([ -1 1 ])
set(gca,'ytick', -1:0.5:1 , 'FontSize',15)
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{32}(V) and \\Delta_{32}(V)'),'FontSize',15)
saveas(gcf, 'C_32_w_original.jpg')


figure
plot(all_years, C_43_w_Gumbel, 'b', 'LineWidth',3)
hold on
plot(all_years, C_43_w_Gumbel_original, 'k', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
box on
xlim([1940 1967])
ylim([ -1 1 ])
set(gca,'ytick', -1:0.5:1 , 'FontSize',15)
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{43}(V) and \\Delta_{43}(V)'),'FontSize',15)
saveas(gcf, 'C_43_w_original.jpg')




%% Figure 6 in the paper

C_21_w_min=C_21_w_R1(:,1);
C_21_w_min(C_21_w_min==-Inf)=[];
if isempty(C_21_w_min)
    C_21_w_min=min(C_21_w_Gumbel)-1;
else
    C_21_w_min=min(C_21_w_min)-1;
end
C_21_w_max=C_21_w_R1(:,2);
C_21_w_max(C_21_w_max==Inf)=[];
if isempty(C_21_w_max)
    C_21_w_max=max(C_21_w_Gumbel)+1;
else
    C_21_w_max=max(C_21_w_max)+1;
end


%%
C_32_w_min=C_32_w_R1(:,1);
C_32_w_min(C_32_w_min==-Inf)=[];
if isempty(C_32_w_min)
    C_32_w_min=min(C_32_w_Gumbel)-1;
else
    C_32_w_min=min(C_32_w_min)-1;
end
C_32_w_max=C_32_w_R1(:,2);
C_32_w_max(C_32_w_max==Inf)=[];
if isempty(C_32_w_max)
    C_32_w_max=max(C_32_w_Gumbel)+1;
else
    C_32_w_max=max(C_32_w_max)+1;
end


%%
C_43_w_min=C_43_w_R1(:,1);
C_43_w_min(C_43_w_min==-Inf)=[];
if isempty(C_43_w_min)
    C_43_w_min=min(C_43_w_Gumbel)-1;
else
    C_43_w_min=min(C_43_w_min)-1;
end
C_43_w_max=C_43_w_R1(:,2);
C_43_w_max(C_43_w_max==Inf)=[];
if isempty(C_43_w_max)
    C_43_w_max=max(C_43_w_Gumbel)+1;
else
    C_43_w_max=max(C_43_w_max)+1;
end


y_axis_min=-6;
y_axis_max=6; %as for men

%% 21

for h=1:n_years
    if  C_21_w_R1(h,1)==-Inf
        C_21_w_R1(h,1)=y_axis_min;
    end
    if  C_21_w_R2(h,1)==-Inf
        C_21_w_R2(h,1)=y_axis_min;
    end
end

for h=1:n_years
    if  C_21_w_R1(h,2)==Inf
        C_21_w_R1(h,2)=y_axis_max;
    end
    if  C_21_w_R2(h,2)==Inf
        C_21_w_R2(h,2)=y_axis_max;
    end
end




figure
plot(all_years, C_21_w_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=C_21_w_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_21_w_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=C_21_w_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_21_w_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, C_21_w_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
box on
xlim([1940 1967])
ylim([ y_axis_min y_axis_max])
set(gca,'ytick',y_axis_min:3:y_axis_max, 'FontSize',15)
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{21}(V)'),'FontSize',15)
saveas(gcf, 'C_21_w.jpg')


%% 32

for h=1:n_years
    if  C_32_w_R1(h,1)==-Inf
        C_32_w_R1(h,1)=y_axis_min;
    end
    if  C_32_w_R2(h,1)==-Inf
        C_32_w_R2(h,1)=y_axis_min;
    end
end

for h=1:n_years
    if  C_32_w_R1(h,2)==Inf
        C_32_w_R1(h,2)=y_axis_max;
    end
    if  C_32_w_R2(h,2)==Inf
        C_32_w_R2(h,2)=y_axis_max;
    end
end


figure
plot(all_years, C_32_w_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=C_32_w_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_32_w_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=C_32_w_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_32_w_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, C_32_w_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
box on
xlim([1940 1967])
ylim([ y_axis_min y_axis_max])
set(gca,'ytick',y_axis_min:3:y_axis_max, 'FontSize',15)
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{32}(V)'),'FontSize',15)
saveas(gcf, 'C_32_w.jpg')


%% 43

for h=1:n_years
    if  C_43_w_R1(h,1)==-Inf
        C_43_w_R1(h,1)=y_axis_min;
    end
    if  C_43_w_R2(h,1)==-Inf
        C_43_w_R2(h,1)=y_axis_min;
    end
end


for h=1:n_years
    if  C_43_w_R1(h,2)==Inf
        C_43_w_R1(h,2)=y_axis_max;
    end
    if  C_43_w_R2(h,2)==Inf
        C_43_w_R2(h,2)=y_axis_max;
    end
end


figure
plot(all_years, C_43_w_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline= C_43_w_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_43_w_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=C_43_w_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_43_w_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, C_43_w_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
box on
xlim([1940 1967])
ylim([ y_axis_min y_axis_max])
set(gca,'ytick',y_axis_min:3:y_axis_max, 'FontSize',15)
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{43}(V)'),'FontSize',15)
saveas(gcf, 'C_43_w.jpg')


